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SUMllARY 


This report is a compilation of the research in orbit determination, 
navigation and optimization in space trajectori es carried out by 
Analytical Mechanics Associates, Inc. under contract to the Lyndon B. 
Johnson Space Center, covering the period October 1, 1975 through 


September 30, 1979. 


J 


TABLE OF CONTENTS 


SUMMARY iil 

INTRODUCTION 1 

I. INITIAL CARTESIAN COORDINATES FOR RAPID PRECISION 

ORBIT PREDICTION I-l 

n. CURVILINEAR PROJECTION DEVELOPMENTS II-l 

III. PERTURBATION-MAGNITUDE CONTROL FOR DIFFERENCE- 

QUOTIENT ESTIMATION OF DERIVATIVES III-l 

IV. DETERMINATION OF THE ACCELEROMETER BIAS FOR 

IN-ORBIT SHUTTLE TRAJECTORIES IV-1 


V 


INTRODUCTION 


Analytical Mechanics Associates, Inc. , under contract to the Lyndon B. Johnson 
Space Center, acted in the capacity of consultants in the areas of orbit determination, 
navigation, optimization techniques and trajectory design for manned space flights. 

In this capacity, several reports were generated and are included in the text of this 
final report. 

(1) Initial Cartesian Coordinates for Rapid Precision Orbit Prediction 

This report contains the equations of motion for a variation of parameters precision 
trajectory prediction which is free of round-off error over long time periods and 
permits large time steps for rapid computation. 

(2) Curvilinear Projection Developments 

This report reviews various methods for accelerating convergence in optimization 
methods using search routines by applying curvilinear projection ideas. 

(3) Perturbation- Magnitude Control For Difference-Quotiert Estimation of Derivatives 
This report develops estimates for choosing perturbation step-sizes used in difference- 
quotient approximations of derivatives for optimization programs. 

(4) Determination of the Accelerometer Bias For In-Orbit Shuttle Trajectories 
This report develops a closed form solution for estimating accelerometer bias when 
using "delta V” accelerometer output for on-board computation of position and 
velocity. 
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SUMMARY 


This report developes the Initial cartesian coordinates of the classical two 
body problem as perturbation parameters for use In rapid and precise computation 
of satellite coordinates in the presence of complex perturbations. The method is 
competitive with the KS theory in that it is numerically stable over long integration 
periods, permitting large integrator step sizes (of the order of 1/50 of the orbital 
period). Moreover, it requires only 7 first order differential equations and re- 
quires no transformation to utilize existing cartesian coordinate perturbation com- 
puting machine routines in the integration of the differential equations of motion. 


I. 


INTRODUCTION 


The numerical Integration of a differential equation suffers from two major 
sources of numerical error. First, an improper choice of the Integration scheme 
may produce spurious solutions in which the numerical error increases exponentially. 
Such errors are due to the characteristic roots of a class of linear differential equa- 
tions with constant coefficients and are associated with the difference equations 
formulae used in the integration scheme. Thus, these errors are independent of 
the actual differential equations whose solution is required and may be controlled, 
or eliminated, by the proper choice of integration scheme and integration step size. 
References [l, 2, and s] contain discussions of these effects. The second error 
source is due to the differential equation itself and is discussed in Reference [l] . 
This error growth comes about from the nature of the solution of the variational 
equations. If the homogeneous variational differential equation contains only 
bounded solutions, the numerical solution of the differential equation will be stable 
and integration can be accurately carried out over long time periods. 

For the differential equations of the satellite numerical instability arises 
from the existence of mixed secular terms in the variational equation of the classical 
two body problem. Reference [4] contains a good illustration of this effect. The 
KS (Reference [s]) eliminates the mixed secular terms in the equations of motion 
by reducing the two body problem to the problem of the linear oscillation. However, 
the computing cost is somewhat increased by requiring a transformation between the 
cartesian space and the regularized variables. Moreover, nine differential equations 
are required in place of the conventional six. The method outlined in this report 
eliminates the mixed secular terms in the variational equations and retains the 
cartesian coordinate frame and time as the independent variable. Only seven 
differential equations are required and some computing time and data storage may 
be saved in comparison with the KS theory for the same accuracy. 
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THE INITLAL CONDITION CARTESLAN ELEMENTS 


The equation of motion of a satellite in the inertial reference frame of the 
earth is given by 


zlL 


R ^ F 


( 1 ) 


dt 


where F represents the pertubation forces other than the central attraction of the 
earth. We seek to represent the instantaneous position and velocity of the satellite 
as an osculating ellipse. Thus, if R(t) and R(t) are the position and velocity we 
wish to describe, the solution is required in the form. 


and 


where 


R(t) = f(0)R^(t) + g(e)R^(t) 


f(0). g(9), ye). gj(e) 


(2a) 

(2b) 


are the classical f and g functions of the two body problem given in References 
[ 6, 7, and S] 

(,e, . 1 . 


g(®) = 


rvasinS da(l-cos6) 

.2 

^/5 " 


, (9) = 

t' ' rr 


(3) 


where 


« a(l-cos0) 

g^,9) . 1 - — ; 


1 = _ 
a r 


^ R * R 
2 0 o 2 


R • R 


r = R 


I -2 


d =R • R 
0 0 0 


(3a) 

coat. 


r = |r I = a(l-cos 9)-^r cos 9-—^ v'^ sin 6 


V/i 


and 6 Ls the difference in eccentric anomaly 

e = E - E 


(3b) 


Given any function b(K^, R^, t) u-e define the total time derivative of h to 


be 
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The first term represents the change in h when the initial conditions are considered 
constant in time, and the sum of the second and third term are che perturbation deri- 
vative of h in which only the variations in the initial condition parameters are per- 
mitted. Thus, 



Bh 

dR 
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9R 
0 
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(4a) 


For R(t), R(t) defined in Eq. (2a) and 2b) to be a solution of £q. (1), we 


require 


and 


^^R(t)=0 

IfBm = F 


(5a) 


fob) 


Equations (5a) and (5b) provide us with the necessary tools to obtain :he 
perturbation derivatives of R^(t) and R^it). 

Since we have 

f,9) {T - g(5) f .9) = 1 (' 

I w 
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Equations (2a) and (2b) may be solved for H^ll) ^ind R^lt) In terms of R(t) 
and R(t). 


Ro(t) = gj(0)R(t)-g(t)R(t) 

(7a) 

Ro(t)*-ft(9)R(t)+f(0)R(t) 

(7b) 


In order to utilize Equations (5a) and (5b) to obtain the perturbation dlfferen- 

tlal equations for the vectors R (t) and R (t) It is necessary to express the f and 

o o 

g functions in terms of the vectors R(t) and R(t). We have, 

i= R • R 
a r /i 

n. o R • R (sin 6) 

r =a(l - cos 6) + r cos 0- _ ‘ (8) 

V/i 

V'M H 

To eliminate the mixed secular terms In the variational equations we follow 
the recommendation of Reference [9] and set the perturbation derivative of the 
change in eccentric anomaly, 6, to zero. Thus, 


^ 9=0 ( 9 ) 

This ensures that the perturbation differential equation will contain period terms and 
not mixed secular terms. Differentiating Equations (7a) and (7b), we have 


Eliminating R(t) and R(ti in Equations ilOai and ilObt we have, 
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Since all the functions are given in terms of 0 and not the time, t, it is 
necessary to have the value of 0 as a function of the independent vr riable t. 
Let 0 be given by 


0-6 V a 


Then, the total derivative of 0 with respect to time is given by. 
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Since 


we have, 
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dt r 
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Integrating Eq. (16), we obtain /3(t). The required 0(t) is given by, 

B (t) 


0(t) = 


Va(t) 


(17) 


The seven differential equations that are to be integrated are Equations (11a), 
(lib) and (16). 


• • 

To obtain the perturbation vector F(t, R, R), we compute R(t) and n(t) 
from Equations (2a) and (2b) (given R^(U. ^(^)) and employ existing 

perturbation routines for F as a function of t, R(t), and R(t). 

III. .A COM.MENT 


It is to be noted that the equations outlined in this report are valid only for 
satellites. Thus, for parabolic and hyperbolic orbits, a univers.al formulation of 
these equations are required which are outlined in Reference [4] . 
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CURVILINEAR PROJECTION DEVELOPMENTS* 
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ABSTRACT 

Gradient Projection is a powerful algorithm for minimization of a function 
subject to constraints (Refs. 1-5), at its best when the constraint functions are 
linear or nearly so. Constraint nonlinearities hamper projection computations, 
often requiring termination of a one-dimensional search in the projected negative 
gradient direction short of the 1-D minimum sought, on account of build-up of 
constraint violations. The constraints must then be restored before another pro- 
jection cycle, at a certain computational expense. The restoration steps taken in 
ihe process of following nonlinear constraint surfaces can be used as a guide to the 
construction of a curve which more nearly follows the constraints than does the 
straight line in the projected gradient direction. This scheme, termed "curvilinear 
projection", was explored in Refs. 7 and 8. The study presently reported carries 
out some computational experiments using a related version of the technique. 

Some other details of projection computations which turn out to be practically im- 
portant are taken up: rules for updating the variable metric in projection when 
early termination of the 1-D search on constraint violation occurs; active-constraint 
logic for screening inequalities that makes use of the Kuhn- Tucker necessary con- 
ditions. Computational comparisons on simple problems are presented. 
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INTRODUCTION 

Original and variable-metric versions of gradient projection algoritIiIns^o^ 
constrained minimization of a function are reported in Refs. 1-3. The present pa- 
per presents some recent improvements, and further investigations of a curved- 
search feature explored in Refs. 7 and 3 which affords improved constraint 
following. 

The resurgence of interest Ln projection, on the part of the present writers, 
came with a surprise in the results of a comparison involving a seemingly slight 
modification of the Kelley-Speyer projection algorithm (Ref. 3). The modification 
was a provision for early updating of the variable metric whenever a screening 
test is passed. A notable convergence Improvement was realized, resulting in the 
projection algorithm, which had been carried along merely for comparison, out- 
performing a more complex algorithm utilizing linear and quadratic penalties. 

The algorithm will first be reviewed in its original equality-constraint ver- 
sion, then the updating rule just mentioned taken up. The restoration of constraints 
and the handling of inequality constraints will be discussed. Attention will then 
turn to the use of search along a curve, proposed in Refs. 7 and 9 with the idea of 
staying closer to constraint surfaces. Some computational e.xperiments will then 
be described. 
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VAP.L\BL E-METRIC PROJECTION 


The projection version of the Davldon-Fletcher-Powell algorithm (Ref. 6) de- 
scribed in the following is essentially the algorithm of Ref. 3; some details are dl^erent, 
however, and the differences Important computationally. The process begins with con- 
straint restoration, usually requiring several cycles; then optimization cycles alternate 
with restorations, which sometimes require more than one cycle. The present section 
will deal with optimization cycles, the following one with restorations. 




The update is performed ooiy if 



Ax^(Af + Ag„X)>0 (5) 

which assures positive definiteness of the updated H. This represents a departure 
from earlier versions of the algorithm (Refs. 3,5) in which termination of one- 
dimensional search on a minimum of f + $ X was required before updating of H 
was permitted, i. e. , updating was deferred until the vicinity of the constrained 
minimum had been reached. 




CONSTRAINT RESTOR.\TION PHASE 


The Initial nulling out of constraint functions often proves more challenging 
than subsequent restorations in that the constraint violations to be dealt with are 
ordinarily larger in magnitude. For this purpose, minimization of a function f 
is employed: 



j=l 


f)-h(f-f^) 


(6) 


This is a weighted sum of squares of the constraint functions plus a term intended 
to counter gross increases in f. The term corresponds to penalty -function treat- 
ment of an inequality ^ - f 5 o . Here h is the heaviside unit step function. The 

k. are determined from 
J 


k. 
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m 



m 





- ,m 
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where k Is Input. This choice would make equal the contribution of each eqxiality 
constraint to the second directional derivative of (6) in its own gradient direction 
S. = 0 , if the constraints were linear. The constraint f - f ^ 0 is Included 

j o 

quadratic-penalty>wise in (6) only during the first restoration sequence, with a 

coefficient k taken as 1/10 the smallest of the k, calculated from (7). The 
o j 

constant ^ is estimated as the Initial value of f g X. 


The metric employed in correction sequences may be denoted A (to dis> 
tinguish it from H of the optimization cycles). It is adjusted approximately for 
changes in the k. , one at a time, using 


A + AA = A - 


Ak. 


1 + Ak g A g 


Ag g/A 

X X 


( 8 ) 


This correction, from Ref. 9, is based on the idea that A approximates f 

XX 

The metric to start the first correction sequence is obtained as A AA from (3), 
using A-i and Ak^ = k^ - 1 (k^ from (7)). If n or more updates are completed 
in this sequence, the emerging DFP metric is carried over to the next: if not, the 
initial metric is carried over. In either case, adjustments for any changes in the 
k^ are performed via (8) before use. Negative increments Ak^ are limited in 
magnitude to insure that the denominator of the fraction in parenthesis does not 
nearly vanish. 


The second and subsequent restoration sequences employ 


T -1 

A.X = - ^®x ®x^ = 


(9) 


together with a one-dimensional search versus a for a minimum of f given by 
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(6), but with the last term deleted. This correction scheme, with ct « 1 and 
without a search, was originally proposed by Rosen (Ref. 1) ; It effects restora- 
tion In a single step for linear g. The existence of the Inverse In (9) (and In (3)) 
requires that the matrix g^ have rank m . This condition Is met at the constrained 
minimum In the classical normal case In which the tangent -plane approximations to 
the constraints are well-defined and distinct. Note that there Is no guarantee that 

A 

(9) Is a direction of descent for f, with general k. values; thus the one-dlmenslonal 

J A 

search may fall and reversion to DFP minimization of f become necessary. 


The magnitude of constraint violation upon which optimization cycles are 
terminated short of a one-dimensional minimum is g^ , where g^ is a pre- 
conceived tolerance and c^ , usually ^ 1 , is a factor adjusted with the aim of 
just permitting restoration with a single cycle of (9), to within the tolerance. Use 
of a single c-factor for all constraints met with only limited success, so a c- 


vector was resorted to, the components adjusted adaptively if somewhat heuristl- 
cally in the following way: c. is increased 10% if a single restoration proves 
successful; it Is halved If two restoration cycles are required; and it Is cut to 


one-quarter If there are additional cycles. 


TREAT^IENT OF TNEQUALITTES 

It is of interest to determine a minimum subject to a mix of equality and 
inequality constraints, the latter expressed by 

g^ a 0 j = m-1, , m-p (10) 

During the initial correcxlon sequence, these are handled penalty- function fashion 
(Ref. 2) , the function i to be minimized given by 


II-<) 
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* i t V] 
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-r ( f - f f 
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h(f-f^) 


( 11 ) 


with the k. determined as though all constraints were equalities: 
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k. 




m+p 

z 

l»l 


I 

f I 


i, 


j = l,2,--,m+p 


(12) 


The determination of the active constraint set for optimization and restora- 
tion cycles proceeds first by excluding those satisfied with a margin g. ^ g . , where 
g^ > 0 is a preset threshhold. Those candidate inequality constraints for which 
gj < Sj » screened further via the Kuhn-Tucker conditions X. s o (Refs. 10 

and 11), using (3) first with ^ the candidates Included, then successively with Kuhn- 
Tucker violators dropped, as many times as necessary, until all s 0 or ail 
candidates are screened out. Inactive constraints are treated in penalty-function 
approximation. 


The Kuhn-Tucker conditions employed apply to the problem of minimizing a 
linear approximation to the function f subject to linearized constraints and to a 
quadratic constraint on step size. They become identical to the Kuhn-Tucker con- 
ditions for the original problem when evaluated at the constrained minimum sought. 


The Kuhn-Tucker screening has generally been found to be worth the com- 
putational expense in reducing tendencies of constraints to switch between active 
and Inactive status from cycle to cycle. The present effort has proceeded on the 
assumption that vector-matrix operations are cheap computationally in relation to 
the cost of gradient and function samples; this is realistic for the trajectory op- 
timization applications of particular interest to the writers. 
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CURVED SEARCH 


Constraint nonlLnearltles hamper projection computations a great deal In 
applications work, often requiring termination of a one-dimensional search short 
of the 1-D minimum sought on account of constraint -violation build-up. It Is of 
Interest to deflect the search away from the straight line In the negative projected 
gradient direction so as to follow approximately the nonlinear constraint Intersection, 
as proposed in Refs. 7 and 8. An improved version of the curved-search technique 
Is given In the following. 

It Is assumed that at least one projection cycle has already been completed 
(the first is done with a linear search) and that the derivative of f g X with re- 
spect to the step-size parameter a has been reduced in magnitude by no more than 
half, that the constraints have been restored by one or more correction cycles, and 
that there has been no change In the active constraint set. 

A curvilinear-projection cycle proceeds by 


Ax = + 7)a 

(13) 

which replaces (11). Here 


C - 

(14) 

Ax - ? 5 
^ “ - 2 

(15) 
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y A.X. 

L -I i 
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a 5 (IG^ 

<y 
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The vector Ax Is the difference between x from the beginning of the preceding 
projection cycle to the present restored point, the beginning of the next. The scalar 
a < 0 is such that the earlier point is regenerated when a = a is introduced into 
(12). Thus (13) generates a parabola in x space which passes through both restored 
points and is tangent to the projected gradient vector at the later one. 

The curved-search sequencing presently in use provides for a possible 
curved-search on all optimization cycles except the first, which uses a linear 
search. Subsequent optimization cycles use a curved search provided the H-update 
test (5) was met on the preceding cycle, none of the inequality constraints has 
changed status (from or to active), and the preceding one-dimensional search did 
not proceed more than halfway to a minimum, as measured in terms of the mag- 
nitude of the derivative of f + g\ with respect to the step-size parameter a . 
Earlier exploratory computations were more cautious in the use of curved searches, 
and generally less successful. The curving steps do nothing beneficial for con- 
jugacy In the subspace of the constraint intersection, but this is already a lost 
cause with OFF when full steps to 1-D minima are not being taken. 


TEST PROBLEMS 


The test problems employed for experimentation were: 


f 


^2 




X. 
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The coefficients were 


a = -10~^ , a. » 10' 


» 1 , » 10^ , b^ = lO'^ 


c » 10 


rl 


The starting point for the numerical computations to be presented was 

« 10 , X 2 » 5 , Xg » 10 

Test Problem had a single equality, = 0 ; ?2 hi; J two equalities, 
® 0 ! So “ ® ■’ equality, g^ = 0 , and one inequality g,, s 0 ; #4 two 

inequalities, g^ ^ 0 , g,, ^ 0 ; #5 two ’’nequalities with the second reversed, 

S, ^ 0 « “So equality, g. = 0 , and one inequality, also reversed. 


COMPUTATIONAL COMPARISONS 

The following results illustrate various features in the context of equality 
constraints. 


NUMBEH OF GR.\DIENT AND FUNCTION-SA^.IPLE COMPUTATIONS 
REQUIRED FOR VARIOUS TEST PROBLEMS 

Test Problem 



]' Original Kelley-Speyer (Ref. 3) 

I 

I Kelley-Speyer plus early H update 


44 ii 691 

j 72 i 1022 

3c ii 4c0 

1 40 & 499 



O 1 




I 

1 


i 


Improved restoration logic added 
Curved search added 


1 X 0*1 Q 


t 


24 i 



■ i 


Linear versus curvilinear projection results are as follows for a larger 
set of test problems: 


NUMBEK OF GRADIENT AND FUNCTION-SAilPLE CCMPUTAHONS 
REQUIRED FOR VARIOUS TEST PROBLEMS 


Test Problem 

41 

#2 

#3 

#4 

45 

wm 

Kelley-Speyer 

29&33S 

31&224 

27&304 

27&310 

274324 

274329 

improved, linear 







Kelley-Speyer 

17&218 

24&2S3 

23&2S9 

24&2S4 

264296 

254235 

curvilinear 


1 


1 




It is noted that the improvement provided by the curved-search feature is considerably 
smaller in problems which Include Inequalities. 


CONCLUDING RE^L■\RKS 


Several developments and refinements of variable-metric projection have 
been presented including a cur.' ed -search technique for nonlinear -constraint- 
surface following, Improved means for control and correction of constraint vio- 
lations, and screening criteria for active-constraint logic fo*r .i ie with inequalities. 
Projection schemes appear quite promising and worth further development and 
evaluation effort. Experience with a larger variety of problems applying the 
various features described is of interest. 
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SUMMARY 

A process for adjusting perturbation magnitude for accurate difference- 
quotient estimation of derivatives is described in the following. The process 
is intended to be carried out sequentially, alternating with iterations of a param- 
eter-optimization algorithm. A more complex and computationally -expensive 
scheme for occasional auxiliary use is also described. 
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I 

I 


I- 

i 

1 
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INTRODUCTION 

It has been recognized for some time that accuracy of partial derivatives 
is important for convergence of variable-metric optimization processes (Refs. 1 
and 2). An adjustment scheme, based upon central differences and truncation 
error, is described in Ref. 3. The adjustment process of the next section is 
similar in concept but focuses on agreement between forward and backward dif- 
ference quotients. 

The perturbation-control logic is being introduced into PEACE, NASA- 
JSC's trajectory-shaping computer program currently in use for space-shuttle 
applications. 
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A TECHNIQUE FOR THE SEQUENTIAL ADJUSTMENT OF PERTURBATION 
MAGNITUDE 


First- and second-derivative estimates f and f are given by the usual 

a aa 

central-difference formulas, 

.f =. t i-f - 2f 

a ” 2 60C ’ aa~ ^^2 

Here f * f(a) Is a fimctlon of a scalar parameter a , f(a 6ct) , and 

f*a f(a - 6 a) . If one requires that the magnitude of the difference between 

forward- and backward-difference estimates of f be at most 

a 


a 





il Mg » 

■d Mg < ^ 


Then an analysis accoimtlng for terms through second order In 6 a leads to 
A a a/b 


where 


b 3 2 1 f I if 

‘ aa' 



if 


2 If 2 
‘ aa 


2 If 1 < 

aa' 


b 

b 


L 

L 


for the largest magnitude perturbation which will hold the truncation (nonlinearity) 
error to the specified level. 
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Bounds are imposed upon the candidate perturbation magnitude, 

5a* 3 if 

3 A if A^ i A s A^ 

3 A^ if A^< A 

To relieve any tendency toward thrashing from cycle to cycle, the value used to 
generate the next derivative estimate is the geometric mean of the old and the 
new: 


Sa"*^ = y 5a ‘ 6a* 

The idea of the scheme is to obtain equality of forward and backward dif- 
ferences to a specified number of significant figures, e. g. , an £ of 10 
corresponds to four-figure agreement. The bounds a^ , b^ , A^^ , A^^ are 
obvious safety devices. 

The adjustment process described generally works well if the function f 
is smooth and nonlinearity the dominant source of error, i. e. , random errors 
are relatively small. The process should survive brief encounters even with 
such errors as jump discontinuities in f. However, if some algorithm is re- 
lentlessly driving a toward the site of a modeling weakness, such as one char- 
acterized by a jump in f or in its first derivative, any perturbation-control 
scheme will be hard pressed to cope. 
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The lower bound Aj^ should be chosen to insure against appreciable 
round-off error, but the choice is not an easy one to make a priori. A lower- 
bound adjustment process, based on the number of significant figures agreement 
between reference and perturbed values of f , will next be described. 

The functions R and R measure the agreement between reference 

-7 

and perturbed values of f (e. g. , R = 10 indicates seven-figure agreement): 


+ 

R 3 


f^- f ! 


max ( t f > I f I ) 


R 2 


I f - f I 

max ( 1 f 1 , 1 f 1 ) 






R^j^ 2 max (R , R ) 


an index of agreement between samples, £ , initially input, may be employed 

to determine the lower bound A , which is to be adjusted from cycle to cycle 

L 

by a rule such as the following: 




= 10 

R R 
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take 



If £_> R 


R MN 


and R ^ 1. 1 R • 
MX MN ’ 




_R 

2 


If £ ^ R . take 
R MN ’ 


£ = £ 

R R 


The lower bound on perturbation magnitude which would correspond to 


the index £ is 
R 


^1 = 




60! 


where &a is the perturbation magnitude actually used in the determination of 

R and £,* . To avert undue downward fluctuation in the bound A, , however, 
MN R ^ 

the geometric mean of A, and the current A, value 

Li L 




is the updated value for use in the next cycle. The corresponding update for 


"r 


-H- / ^ 

S ^ '' S' s 


The process just described adjusts the index of agreement and the 
perturbation bound A, upward or downward whenever small-perturbation samples 

i-i 

become available. It attributes any disagreement between forward and backward 
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differences to random error, adjusting £_ upward when agreement Is poor, 

xi 

downward when agreement Is good. The factors 1.1, 10 , and . 5 are rather 
arbitrary and intended merely to be suggestive. 


^ 1 4 * 

An alternative scheme for the calculation of ^ and will next be 
described. This is in process of evaluation at the present writing. The round- 
off error magnitude is estimated from the residuals of a least-squares fit In 
the vicinity of the minimum found in the preceding one-dimensional search. 

Of main interest for fitting are the samples in the band k . = 6k , where 


n 



corresponds to the perturbation magnitude employed in the most recent estimate 

of the gradient vector f^ . All of the samples in the band are used for the fit 

provided there are at least four in addition to the k . sample. If not, more 

min 

close samples are added to make a total of five. If any of these fall outside the 
band k^^^ = 10 6k , the computation is abandoned and no update of and 
Ay carried out. 

After the least-squares fit has been completed, the residuals at k^^^ 

and the two closest points are examined. If these all have the same sign, the 

attempt to update A^ and aT^ is abandoned, otherwise the average of the 

absolute values is adopted for £_ . Note that £ in this mode is a scalar, 

K R 

applicable to adjustment of all components. The component -by-component up- 

_9 

date proceeds as follows: If f ^ 10 , take 

^ aa 


IIH5 





otherwise do not update A. for the particular component. 

Xj 

the geometric mean: 


Obtain from 





While the parameter a Is a scalar In the preceding formulae, the pro- 
cess Is intended to apply to each component of the parameter vector in turn. 

The idea of tailoring the choice of pertuilDation magnitude to the particular 
component is hardly earth-shaking but, in fact, it is not often carried out in 
practice. 
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A SEARCH TECHNIQUE FOR PERTURBATION-MAGNITUDE DETERMINATION 


Applications of variable-metric optimization algorithms to complex 
models sometimes encounter convergence difficulties which are attributed to 
numerical errors, real or fancied. The analyst, In a moment of paranoia, 
suspects overlap between the regions of round-off and truncation, with no good 
compromise choice of perturbation magnitude for the generation of secant par- 
tlals available. The situation may be thought serious enough to warrant a 
search, usually a tedious cut-and-try affair. The present section describes 
a mechanization of such a search, employing the one-dimensional minimiza- 
tion technique of Ref. 4. 

A function Q may be defined which measures the error between for- 
ward and backward derivatives: 



Here the denominator-normalization choice between fonvard and backward 

derivative magnitudes is intended to avoid small -divisor difficulties. The 

2 

proposed search is for a minimum of Q subject to a round-off constraint 
on perturbation magnitude introduced via quadratic penalty. Thus 

mm [Q + h(£^^- , 

where h , the Heaviside unit step function, is unity for argument ^ 0 , zero 
otherwise. The function R,,., is that of the previous section. The round-off 
constraint incorporated via penalty is R^^^- - < 


A choice of agreement-index Is no easier to make than in the case 
of the sequential-adjustment process. Values obtained in the course of sequential- 
adjustment cycles may be employed, or, in difficult cases, a range of values 
investigated. 
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CONCLUDC^G REMARKS 


The two processes described herein have been given only limited trials 
and require further attention. An initial application of the adjustment scheme 
indicates, a bit surprisingly, that it may turn out to have particular merit as 
a "debugging" tool. The first hill -scale application presently planned is to the 
accelerated-gradient trajectory shaping program PEACE at NASA-JSC. 
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SUMMARY 


This report yields a closed form deterministic solution of the accelerometer 
bias for in-orbit shuttle trajectories given the difference between the ground 
tracking solution of the trajectory and the on board esUniate of the trajectory 
using platform accelerometers. 


INTRODUCTION 

Due to the presence of small unpredictable accelerations acting on the 
shuttle in orbit it may be necessary to use the stable platform accelerometers 
to integrate the vehicle state. In this mode the most significant error will be 
the accelerometer bias. A method is provided whereby these biases may be 
estimated by using the difference between the vehicle inertial position vectors 
as measured by ground tracking and that obtained by numerical integration of 
the on board accelerometers. 
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SOLUTION OF THE ACCELEROMETER BIAS IN ORBIT 


Let the inertial equations of motion of the orbiting vehicle as measured by 
the ground tracking system be 




( 1 ) 


Let the inertial equations of motion of the orbiting vehicle as measured by 
the on board accelerometers be given by 


' * ^2 

= -M — - - B 


( 2 ) 


where is the true specific force acting and B an inertial vector of the small 
accelerometer bias error. The difference between the inertial vector positions 
is assumed small compared to either or . 

That is. 


and 




Rj-Rj 




R^l >> .001 


R 1 > > .001 


( 3 ) 


Under these conditions the difference in the orbiting states is given by 


-Rj-Ri 


«2 - \ 


= o(t.y 


\ - ''I'l 


‘ 10 

+ \ o (t, T) ( ) dr 




B 


( 4 ) 


The constant accelerometer bias vector \ B / is given by 

"o Ri(t)> 


B = (C) 




t 


B R^ (t) 


\ (toV 




j 


(5) 
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The 3 matrices ( C ) 


\hR(iy VaR^ito)/ 


are given by 


(C) = R^ (t) R^'^(t) + R^(t) R^{t) + Cg R^(t) R^'^(t) 


1 T 


+ R^ (t) R^ (t) + Cg I (3) 


(6) 


a Ri(t)' 


Hi (to)' 


1= R^ (to) R^'^(to) + R^(to) r \) - 83 R^(t^) R^^(to) 


+ R^ (to) R (to) + a^ I (3) 


a R^(t)' 


1 = ’'l ^ <‘o> "'=2 ^<‘c> *'=3 ^1 “o' 


'T 


+ R^ (to) R^(to) - bg I (3) 


The matrix I (3) is the 3x3 unit matrix 


1(3) = 



(7) 


The scalar coefficients a. , b. and c. are given by 
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= — 



(2 ^ G 3 ) 


^2 = 


S 


r r V U 

0 


( 3 ) 
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(8 con' t. ) 



where 


2C.-$0^* 'o ' °3 ■ ^ °2 ’ 



(2G4-SG3' 


= lRi(‘o>l 

r = I (t) j 



Rl(to) 

(t) ‘ 

R^^ (t) 

R^ (to) 

• Ri (to) 

R^ (t) • 

Rl (t) 

/m (t - 

to) 


— » -r — 



(3a) 


IV - 3 


G. = 


^2 = 




(Y ^ — 

(2k+ 1) 1 


G G 
12 

vVoS 


Hi . 


M rro 
2 


M^o 


r^ 


"4 = 


b. = 

D 


SPi 

372 

rji 


§72 


r G- 

0 1 


d G- 
o 2 


y\) hj - ( d Y h^ + h.) h)/ Ah. 

(< S ' '’4'» hj - ( hj * hj d + h.) h^7 Ah. 


d h^) - ( d + h j ) / Ah^ 


3' 4 


'4 5' 2 


“ 0 '’2 * *>4 <*> >>2 - ( hj + hj d * h.) h^^/ Ah. 


c. = 
o 


1/h. 

D 


where 
h = 


2 r y ‘ 


(3a con*t* ) 


(Sb) 


(8c) 


(Sd) 
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K = — (-2G + 


°o <=5 ^ . , 

+ + G ) 

2 2 12 


"3 = 


S <2G4-^G3)-2 G3G^ 




3/21 


(8d) 


1 . 3 ^ .1 ^ 5 , 


^.3(2G^-3/2^Gg.-G^ Gg . ~ ^ S " 48 ^ S> 


T/2 (°2 <°3 - ^°2> ^ °2°3 * 2 G. - 1 - 4-G^ - -fjG^ 


.) 


h. = 


“ I 3 (2 G, - 3/2 8 G, * i G^ G. 4 G^ . -^ b \) 

--T/2 (SP°5-^°4>-S<2°4-^°3>) 


h, = 


r . G 
M 


/ ^2 2 N 

^2 2 ^ 2 Vm ^ ° 


/* 1 R ^ 

- 2 G +— G G +— G + -£- G 
5 2 0 5 2 4 12 2 


~2 ^ 


+ r (G^G2 - 1 G^ G3 - I SG^) 




^ ■— ^^2 

2 */M 




A = ( h3) ( r" v^ - d“) + h. (h^ r“ + h3 d + h^ d + h.) 


Since the vector bias, { B ) , is considered to be constant over an estimate 

^ ^ f 1 

period, it is recommended that several estimates of < B > be obtained using the 
same initial epoch ( tQ , (Iq) , (t^) , R2 (t^), R2 (to) ) nnd several different 
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later times, t 




. The recommended value of 



will be given by 



( 9 ) 


j 


IV - 6 


